﻿using System;
using System.Collections.Generic;
using System.Text;
using SocialExplorer.IO.FastDBF;
using System.IO;

namespace Soil
{
    class SoilLayer
    {
        private int _layer_no;
        private double _sol_z;
        private double _sol_bd;
        private double _sol_awc;
        private double _sol_k;
        private double _sol_cbn;
        private double _clay;
        private double _silt;
        private double _sand;
        private double _rock;
        private double _sol_alb;
        private double _usle_k;
        private double _sol_ec;
        private double _sol_cal;
        private double _sol_ph;

        private double _cec;
        private string _sol_id;
        private double _sol_top_z;

        public string SOL_ID { get { return _sol_id; } }
        public int Layer { get { return _layer_no; } set { _layer_no = value; } }
        public double SOL_Z { get { return _sol_z; } set { _sol_z = value; } }
        public double SOL_BD { get { return _sol_bd; } set { _sol_bd = value; } }
        public double SOL_CBN { get { return _sol_cbn; } set { _sol_cbn = value; } }
        public double CLAY { get { return _clay; } set { _clay = value; } }
        public double SILT { get { return _silt; } set { _silt = value; } }
        public double SAND { get { return _sand; } set { _sand = value; } }
        public double ROCK { get { return _rock; } set { _rock = value; } }
        public double SOL_ALB { get { return _sol_alb; } set { _sol_alb = value; } }
        public double USLE_K { get { return _usle_k; } set { _usle_k = value; } }
        public double SOL_EC { get { return _sol_ec; } set { _sol_ec = value; } }
        public double SOL_TOP_Z { get { return _sol_top_z; } set { _sol_top_z = value; } }
        public double SOL_K { get { return _sol_k; } set { _sol_k = value; } }
        public double SOL_CAL { get { return _sol_cal; } set { _sol_cal = value; } }
        public double SOL_PH { get { return _sol_ph; } set { _sol_ph = value; } }
        public double SOL_AWC { get { return _sol_awc; } set { _sol_awc = value; } }
        public double Thick { get { return _sol_z - _sol_top_z; } }

        public SoilLayer(DbfRecord layerRecord)
        {
            _layer_no = int.Parse(layerRecord["LAYER_NO"]);
            _sol_id = layerRecord["SOIL_ID"];

            _sol_top_z = readDouble(layerRecord, "UDEPTH") * 10;
            _sol_z = readDouble(layerRecord, "LDEPTH") * 10;

            _sol_bd = readDouble(layerRecord, "BD");

            //calculate awc
            double kp33 = readDouble(layerRecord, "KP33");
            double kp1500 = readDouble(layerRecord, "KP1500");
            _sol_awc = SWAT_EMPTY_VALUE;
            if(kp33 > 0 && kp1500 > 0)
                _sol_awc = (kp33 - kp1500) / 100.0;

            _sol_k = readDouble(layerRecord, "KSAT") * 10; //cm/h->mm/h
            _sol_cbn = readDouble(layerRecord, "ORGCARB");
            _clay = readDouble(layerRecord, "TCLAY");
            _silt = readDouble(layerRecord, "TSILT");
            _sand = readDouble(layerRecord, "TSAND");

            //convert volume percent to weight percent            
            _rock = SWAT_EMPTY_VALUE;
            double cofrag = readDouble(layerRecord, "COFRAG");
            if (cofrag > 0 && cofrag <= 100 && _sol_bd > 0)
                _rock = cofrag * 1.4 * _sol_bd;
            if (_rock < 0 || _rock > 100) _rock = SWAT_EMPTY_VALUE;

            //calculate albedo using equation 
            _sol_alb = 0.1; //minimum value for albedo from SWAT code
            if(_sol_cbn > 0 && _sol_cbn <= 1)
                _sol_alb = 0.4 * Math.Exp(0.4 * 1.72 * _sol_cbn);

            //calculate USLE_K
            calculateUSLEK();

            _sol_ec = readDouble(layerRecord, "EC");
            _sol_cal = readDouble(layerRecord, "CACO3");
            _sol_ph = readDouble(layerRecord, "PH2");
            _cec = readDouble(layerRecord, "CEC");
        }

        public static SoilLayer WaterSoilLayer
        {
            get
            {
                SoilLayer layer = new SoilLayer(0, 25.4);
                layer.Layer = 1;
                layer.SOL_BD = 1.72;
                layer.SOL_K = 260;
                layer.SOL_ALB = 0.23;
                return layer;
            }
        }

        public SoilLayer(double top_z, double bottom_z)
        {
            _sol_top_z = top_z;
            _sol_z = bottom_z;
        }

        public void LowerDepth(double adjustDepth)
        {
            _sol_z += adjustDepth;
            _sol_top_z += adjustDepth;
        }

        /// <summary>
        /// Calculate USLE_K based on equation 4:1.1.5~4:1.1.9 in SWAT theory doc.
        /// </summary>
        private void calculateUSLEK()
        {
            if (_sand == SWAT_EMPTY_VALUE || _silt == SWAT_EMPTY_VALUE || 
                _clay == SWAT_EMPTY_VALUE || _sol_cbn == SWAT_EMPTY_VALUE) 
            {_usle_k = 1.0; return; }

            double f_csand = (0.2 + 0.3 * Math.Exp(-0.256 * _sand * (1.0 - _silt / 100.0)));
            double f_cl_si = Math.Pow(_silt / (_clay + _silt), 0.3);
            double f_orgc = 1.0 - 0.25 * _sol_cbn/(_sol_cbn + Math.Exp(3.72 - 2.95*_sol_cbn));
            double sandFraction = 1.0 - _sand / 100;
            double f_hisand = 1.0 - 0.7 * sandFraction/(sandFraction + Math.Exp(-5.51+22.9*sandFraction));
            _usle_k = f_csand * f_cl_si * f_orgc * f_hisand;

            if(_usle_k <= 0 || _usle_k > 1) _usle_k = 1.0;
        }

        public static double readDouble(DbfRecord rec, string col)
        {
            double v = SWAT_EMPTY_VALUE;
            double.TryParse(rec[col], out v);
            if(v == NSDB_EMPTY_VALUE) v = SWAT_EMPTY_VALUE;
            return v;
        }



        public override string ToString()
        {
            return string.Format("{0},{1},{2},{3},{4},{5},{6},{7},{8},{9},{10},{11}",
               _sol_z,_sol_bd,_sol_awc,_sol_k,_sol_cbn,_clay,_silt,_sand,_rock,_sol_alb,_usle_k,_sol_ec);
        }

        public static string EMPTY_VALUES = "0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0";
        public static double NSDB_EMPTY_VALUE = -9.0;
        public static double SWAT_EMPTY_VALUE = 0.0;

        public SoilLayer Clone()
        {
            return this.MemberwiseClone() as SoilLayer;
        }

        public SoilLayer Clone(double depth_up, double depth_bottom)
        {
            if (_sol_top_z >= depth_bottom || _sol_z <= depth_up) return null;

            SoilLayer layer = Clone();
            layer.SOL_TOP_Z = Math.Max(depth_up, SOL_TOP_Z);
            layer.SOL_Z = Math.Min(depth_bottom, SOL_Z);
            return layer;
        }

        /// <summary>
        /// if overlay with given layer
        /// </summary>
        /// <param name="layer"></param>
        /// <returns></returns>
        public bool isOverlay(SoilLayer layer)
        {
            if (SOL_TOP_Z <= layer.SOL_Z ||
                SOL_Z >= layer.SOL_TOP_Z) return false;

            return true;
        }

        /// <summary>
        /// if connected with given layer
        /// </summary>
        /// <param name="layer"></param>
        /// <returns></returns>
        public bool isConnected(SoilLayer layer)
        {
            if (layer == null) return false;

            return this.SOL_TOP_Z == layer.SOL_Z || this.SOL_Z == layer.SOL_TOP_Z;
        }

        public bool isSameDepth(SoilLayer layer)
        {
            if (layer == null) return false;

            return SOL_TOP_Z == layer.SOL_TOP_Z && SOL_Z == layer.SOL_Z;
        }

        /// <summary>
        /// Connect two layers to one. The properties would be calculated based on thick
        /// </summary>
        /// <param name="layer"></param>
        /// <returns></returns>
        public SoilLayer Connect(SoilLayer layer)
        {
            if (!isConnected(layer))  //only connected derivs could be connected
                throw new Exception("Layers are not connected.");

            SoilLayer connectLayer = new SoilLayer(
                Math.Min(this.SOL_TOP_Z, layer.SOL_TOP_Z),
                Math.Max(this.SOL_Z, layer.SOL_Z));

            //calculate height average of all properties
            double ratio_this = Convert.ToDouble(Thick) / Convert.ToDouble(connectLayer.Thick);
            double ratio_another = 1.0 - ratio_this;

            //todo -- add property calculation here
            CalculateProperties(connectLayer,
                this, ratio_this,
                layer, ratio_another);

            return connectLayer;
        }

        /// <summary>
        /// Merge two deriv to one. The properties would be calculated based percentage
        /// </summary>
        /// <param name="deriv"></param>
        /// <returns></returns>
        public SoilLayer Merge(double percentage, SoilLayer layer_other, double percentage_other)
        {
            if (layer_other == null) return this;

            if (!isSameDepth(layer_other))  //only same derivs could be merged
                throw new Exception("Layers are not same size.");

            SoilLayer mergeDeriv = this.Clone();

            //todo -- add property calculation here based on percentage
            CalculateProperties(mergeDeriv,
                this, percentage,
                layer_other, percentage_other);

            return mergeDeriv;
        }

        /// <summary>
        /// Calculate new properties based on two other derivs and ratio
        /// </summary>
        /// <param name="newLayer"></param>
        /// <param name="layer"></param>
        /// <param name="percentage"></param>
        /// <param name="layer_other"></param>
        /// <param name="percentage_other"></param>
        private static void CalculateProperties(SoilLayer newLayer, SoilLayer layer, double percentage, SoilLayer layer_other, double percentage_other)
        {
            newLayer.SOL_BD = layer.SOL_BD * percentage + layer_other.SOL_BD * percentage_other;
            newLayer.SOL_AWC = layer.SOL_AWC * percentage + layer_other.SOL_AWC * percentage_other;
            newLayer.SOL_K = layer.SOL_K * percentage + layer_other.SOL_K * percentage_other;
            newLayer.SOL_CBN = layer.SOL_CBN * percentage + layer_other.SOL_CBN * percentage_other;
            newLayer.CLAY = layer.CLAY * percentage + layer_other.CLAY * percentage_other;
            newLayer.SILT = layer.SILT * percentage + layer_other.SILT * percentage_other;
            newLayer.SAND = layer.SAND * percentage + layer_other.SAND * percentage_other;
            newLayer.ROCK = layer.ROCK * percentage + layer_other.ROCK * percentage_other;
            newLayer.SOL_ALB = layer.SOL_ALB * percentage + layer_other.SOL_ALB * percentage_other;
            newLayer.USLE_K = layer.USLE_K * percentage + layer_other.USLE_K * percentage_other;
            newLayer.SOL_EC = layer.SOL_EC * percentage + layer_other.SOL_EC * percentage_other;
            newLayer.SOL_CAL = layer.SOL_CAL * percentage + layer_other.SOL_CAL * percentage_other;
            newLayer.SOL_PH = layer.SOL_PH * percentage + layer_other.SOL_PH * percentage_other;
        }
    }
}
